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traps by Fisher zeros 



O 

o 

(N 



Oliver Miilken, Peter Borrmann, Jens Harting, and Hcinrich Stamcrjohanns 
Department of Physics, Carl von Ossietzky University Oldenburg, D-26111 Oldenburg, Germany 

(Dated: February 1, 2008) 

We present a detailed description of a classification scheme for phase transitions in finite systems 
based on the distribution of Fisher zeros of the canonical partition function in the complex temper- 
ature plane. We apply this scheme to finite Bose-systems in power law traps within a semi-analytic 
approach with a continuous one-particle density of states Sl(E) ~ E for different values of d and 
to a three dimensional harmonically confined ideal Bose-gas with discrete energy levels. Our results 
indicate that the order of the Bose-Einstein condensation phase transition sensitively depends on 
the confining potential. 
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I. INTRODUCTION 

In 1924 S. Bose and A. Einstein predicted that in a 
system of bosons at temperatures below a certain crit- 
ical temperature T c the single-particle ground state is 
macroscopically occupied |lj]. This effect is commonly 
referred as Bose-Einstein condensation and a large num- 
ber of phenomena, among others the condensation phe- 
nomena in alkali atoms, the superfluidity of 4 He and the 
superconductivity, are identified as signatures of this ef- 
fect. However, the physical situation is very intricate in 
most experiments. 

Recent experiments with dilute gases of alkali atoms 
in magnetic and optical J3| traps are in some sense 
the up to now best experimental approximation of the 
ideal non-interacting Bose-Einstein system in an exter- 
nal power law potential. The achievement of ultra-low 
temperatures by laser cooling and evaporative cooling 
opens the opportunity to study the Bose-Einstein con- 
densation under systematic variation of adjustable ex- 
ternal parameters, e.g the trap geometry, the number of 
trapped atoms, the temperature, and by the choice of 
the alkali atoms the effective interparticle interactions. 
Even in the approximation of non-interacting particles 
the explanation of these experiments requires some care, 
because the number of bosons in these novel traps is finite 
and fixed and the standard grand-canonical treatment is 
not appropriate. The effect of the finite particle numbers 
on the second moments of the distribution function, e.g. 
the specific heat and the fluctuation of the ground state 
occupation number has been addressed in a number of 
publications J| §]. In @, §| we have presented a recursion 
method to calculate the canonical partition function for 
non-interacting bosons and investigated the dependency 
of the thermodynamic properties of the condensate on 
the trap geometry. 

The order of the phase transition in small systems sen- 
sitively depends on finite size effects. Compared to the 
macroscopic system even for as simple systems as the 3- 
dimensional ideal gas the order of the phase transition 
might change for mesoscopic systems where the number 
of particles is finite or for trapped gases with different 
trap geometries. 



In this paper we address the classification of the phase 
transition of a finite number of non-interacting bosons 
in a power law trap with an effective one-particle den- 
sity of states Q,(E) = E d ~ x being formally equivalent to 
a d-dimensional harmonic oscillator or a 2d-dimcnsional 
ideal gas. We use a classification scheme based on the 
distribution of zeros of the canonical partition function 
initially developed by Grossman et al. fiL and Fisher 
et al. ||, which has been extended by us f| as a classi- 
fication scheme for finite systems. On the basis of this 
classification scheme we are able to extract a qualita- 
tive difference between the order of the phase transition 
occuring in Bose-Einstein condensates in 3-dimensional 
traps [llj and in 2-dimensional traps which was re- 
cently discovered by Safonov et al. in a gas of hydrogen 
atoms absorbed on the surface of liquid helium jl2j . Since 
we do not consider particle interactions this difference is 
only due to the difference in the confining potential. 

We give a detailed review of the classification scheme 
in Sec. II. In Sec. Ill we present the method for the calcu- 
lation of the canonical partition function in the complex 
plane and describe details of the numerical implementa- 
tion. Our results for d = 1 — 6 and particl e n umbers 
varying from 10 to 300 are presented in Sec. Ill as well 
as calculations for a 3-dimensional parabolically confined 
Bose-gas. 



II. CLASSIFICATION SCHEME 

In 1952 Yang and Lee have shown that the grand 
canonical partition function can be written as a func- 
tion of its zeros in the complex fugacity plane, which lie 
for systems with hard-core interactions and for the Ising 
model on a unit circle [ p"3[ . 

Grossmann et al. JjJ and Fisher || have extended this 
approach to the canonical ensemble by analytic contin- 
uation of the inverse temperature to the complex plane 
(3 — > B = (3 + it. Within this treatment all phenomeno- 
logically known types of phase transitions in macroscopic 
systems can be identified from the properties of the dis- 
tribution of zeros of the canonical partition function. 

In U we have presented a classification scheme for fi- 
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nite systems which has its macroscopic equivalent in the 
scheme given by Grossmann. As usual the canonical par- 
tition function reads 



Z(B) = J dE Q(E) cxp(-BE), 



(1) 



which we write as a product Z(B) = Zn m (B)Z int (B), 
where Zn m (B) describes the limiting behavior of Z(B) 
for T — > oo imposing that lim-r^oo Z- mt (B) = 1. This 
limiting partition function will only depend on the ex- 
ternal potential applied to the system, whereas Zi nt (B) 
will depend on the specific interaction between the sys- 
tem particles. E.g. for a iV-particle system in a d- 
dimensional harmonic trap Zn m (B) = B~ dN and thus the 
zeros of Z(B) are the same as the zeros of Z- ln t(B). Since 
the partition function is an integral function, the zeros 
Bk = B*_ k — (3k + iT~k (k £ N) are complex conjugated 
and the partition function reads 



Z(B) 



Ziim(B) Z int (0) exp(Bd B \nZ int (0)) 



feeN 



exp 



B B , 



The zeros of Z(B) are the poles of the Helmholtz free 
energy F(B) = —^\nZ(B), i.e. the free energy is ana- 
lytic everywhere in the complex temperature plane ex- 
cept at the zeros of Z(B). 

Different phases are represented by regions of holomor- 
phy which are separated by zeros lying dense on lines in 
the complex temperature plane. In finite systems the ze- 
ros do not squeeze on lines which leads to a more blurred 
separation of different phases. We interpret the zeros as 
boundary posts between two phases. The distribution 
of zeros contains the complete thermodynamic informa- 
tion about the system and all thermodynamic properties 
are derivable from it. Within this picture the interaction 
part of the specific heat is given by 



^ 2 E 



(B k - B) 2 (B* k - Bf 



(3) 



The zeros of the partition function are poles of Cy{B). 
As can be seen from Eq. (||) a zero approaching the real 
axis infinitely close causes a divergence at real tempera- 
ture. The contribution of a zero Bk to the specific heat 
decreases with increasing imaginary part rfc. Thus, the 
thermodynamic properties of a system are governed by 
the zeros of Z close to the real axis. 

The basic idea of the classification scheme for phase 
transitions in small systems presented in Q is that the 
distribution of zeros close to the real axis can approxi- 
mately be described by three parameters, where two of 
them reflect the order of the phase transition and the 
third merely the size of the system. 

We assume that the zeros lie on straight lines (see 
Fig. [I]) with a discrete density of zeros given by 
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FIG. 1: Schematic illustration of the zeros in the complex 
temperature plane. 



with k = 2,3,4, and approximate for small r the 
density of zeros by a simple power law 4>(t) ~ r a . Con- 
sidering only the first three zeros the exponent a can be 
estimated as 



Injurs) - ln</>(T 2 ) 
In r 3 - In r 2 



(5) 



The second parameter to describe the distribution of ze- 
ros is given by 7 = tanz/ ~ — (3\) / '(ra — ti) where v is 
the crossing angle of the line of zeros with the real axis 
(see Fig. [j]). The discreteness of the system is reflected 
in the imaginary part n of the zero closest to the real 
axis. 

In the thermodynamic limit we have always t\ — > 0. 
In this case the parameters a and 7 coincide with those 
defined by Grossmann et al B, who have shown how 
different types of phase transitions can be attributed to 
certain values of a and 7. They claimed that a = 
and 7 = corresponds to a first order phase transition, 
second order transitions correspond to < a < 1 with 
7 = or 7 ^ 0, third order transitions to 1 < a < 2 with 
arbitrary values of 7, and that all higher order phase 
transition correspond to a > 1. For macroscopic systems 
(with n — > 0) a cannot be smaller than zero, because this 
would cause a divergence of the internal energy. However 
in small systems with a finite t± this is possible. 

In our classification scheme we therefore define phase 
transitions in small systems to be of first order for a < 0, 
while second and higher order transitions are defined in 
complete analogy to the Grossmann scheme augmented 
by the third parameter t\. The definition of a critical 
temperature [3c in small systems is crucial and ambigu- 
ous since no thermodynamic properties diverge. Thus, 
different definitions are possible. We define the critical 
temperature as (3 cu t = Pi — yri, i.e. the crossing point 
of the approximated line of zeros with the real tempera- 
ture axis. An alternative definition is the real part of the 
first complex zero /3i. In the thermodynamic limit both 
definitions coincide. 

Comparing the specific heats calculated for different 
discrete distributions of zeros shows the advantages of 
this classification scheme. Fig. || shows (a) three distri- 
butions of zeros lying on straight lines corresponding to 
a first order transition (a = and 7 = 0), a second order 



3 



transition (a = 0.5 and 7 = —0.5), and a third order 
phase transition (a — 1.5 and 7 = —1) and (b) the perti- 
nent specific heats. In all cases the specific heat exhibit 
a hump extending over a finite temperature region and 
cannot be used to classify the phase transition. In con- 
trast, even for very small systems (large ti) the order of 
the phase transition is extractable from the distribution 
of zeros. 




FIG. 2: Plot of (a) generated zeros lying on straight lines 
to simulate first (a = and 7 = 0), second (a = 0.5 and 
7 = —0.5), and third (a = 1.5 and 7 = —1) order phase 
transitions and (b) the appropriate specific heats per particle. 



The zeros of the canonical partition function have a 
distinct geometrical interpretation which explains the 
smoothed curves of the specific heat and other thermo- 
dynamic properties in finite systems. 

Fig. U shows (a) the ground state occupation number 
\r]o(B)\/N in the complex temperature plane and (b) the 
ground state occupation number at real temperatures for 
a finite ideal Bose gas of TV = 120 particles, where i]o(B) 
is given by the derivative of the logarithm of the canonical 
partition function Z(B) with respect to the ground state 
energy e , i.e. 770(B) = -±d t0 Z(B)/Z{B). 




FIG. 3: Comparison of (a) \rjo\/N with (b) the appropriate 
value of 770 at real temperatures for a 120 particle harmoni- 
cally trapped ideal Bose-gas (note that K — /cb = cj — 1). 

Zeros of the partition function are poles of 770(B) and 
are indicated by dark spots, which influence the value 



of the ground state occupation number at real tempera- 
tures impressively. Every pole seems to radiate onto the 
real axis and therefore determines the occupation num- 
ber at real temperatures. This radiation extends over a 
broad temperature range so that the occupation number 
for real temperatures does not show a discontinuity but 
a smoothed curve. A closer look at Eq. (|J) gives the 
mathematical explanation for this effect. The discrete 
distribution of zeros, i.e. t\ > 0, inhibits the specific 
heat and all other thermodynamic properties to show a 
divergency at some critical temperature because the de- 
nominators of the arguments of the sum remain finite. 

Without going into a detailed analysis we note that in 
the thermodynamic limit the parameter a is connected 
to the critical index for the specific heat by 



C v ~ (/? - (3 c ) a ~ 



(6) 



However, since critical indices are used to describe the 
shape of a divergency at the critical point an extension 
to small systems seems to be more or less academical. 

The introduction of complex temperatures might seem 
artificial at first sight but, in fact, the imaginary parts 
Tfc of the complex zeros Bk have an obvious quantum 
mechanical interpretation. We write the quantum me- 
chanical partition function as 

Z((3 + ir/h) = Tr(exp(-?;riJ/ft)exp(-/37J)) (7) 
= ^ can |exp(-ir^'/?i)|* can \ (8) 

= (#can(t = 0)|tf can (i =r)), (9) 

introducing a canonical state as a sum over Boltzmann- 
weighted eigenstates *can) = J2k exp(-/3e fc /2) \4> k ). We 
explicitly write the imaginary part as t/H since the di- 
mension is 1/ [energy] and the imaginary part therefore 
can be interpreted as time. Then the imaginary parts Tk 
of the zeros resemble those times for which the overlap 
of the initial canonical state with the time evoluted state 
vanishes. However, they are not connected to a single 
system but to a whole ensemble of identical systems in a 
heat bath with an initial Boltzmann distribution. 



III. BEC IN POWER LAW TRAPS 

In this section we assume a continuous single particle 
density of states Cl(E) = E as an approximation for 
a d-dimensional harmonic oscillator or a 2d-dimcnsional 
ideal gas. E.g. for the harmonic oscillator this corre- 
sponds to the limit of hw — > and taking only the lead- 
ing term of the degeneracy of the single particle energy 
levels. The one-particle partition function is given by the 
Laplace transformation 



Zi{B) = / dE E d ~ 1 exp(-BE) = (d - 1)! B~ 



(10) 



The canonical partition function for N non-interacting 
bosons can be calculated by the following recursion H 



Z N (B) 



1 N 

- ]T Z 1 {kB) z 



N-k 



(11) 



k=l 



where Z±(kB) = J2i ex P(~^^ e i) is the one-particle par- 
tition function at temperature kB and Z$(B) — 1. For 
small particle numbers this recursion works fine, even 
though its numerical effort grows proportional to TV 2 . 

With ( [To|) as Zx Eq. ([ll]) leads to a polynomial of order 
N in (1/B) d for Zm which can be easily generated using 
Maple or Mathematica. The zeros of this polynomial 
can be found by standard numerical methods. 
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FIG. 4: Distribution of zeros for Bose-Einstein Condensates 
with continuous one-particle density of states £l(E) — E' 1 ^ 1 
for d = 1 - 6. 

Fig. [| displays the zeros of the iV-particle partition 
function for d — 1 — 6 in the complex temperature plane 
for particle numbers N = 25, 50 and 100. For d = 2 — 6 
the zeros approach the positive real axis with increas- 
ing particle number and are shifted to higher tempera- 
tures which is already at first sight an indicator of phase 
transitions. For d = 1 the zeros approach the real axis 
only at negative temperature. This behavior is consistent 
with the usual prediction that there is no Bose-Einstein 
condensation for the one-dimensional harmonic oscillator 
and the two-dimensional ideal Bose gas Efl. 
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FIG. 5: Specific heat scaled by dN of Bose-Einstein Con- 
densates with continuous one-particle density of states for 
d= 1 -6. 



The symmetry of the distributions of zeros is due to the 
fact that Zjy is a polynomial in B~ d . For this reason it 
can be inferred that for d — > oo the zeros lie on a perfect 
circle. 

Fig. U shows the corresponding specific heats calcu- 
lated using equation (||). As expected, for d = 1 the spe- 
cific heat has no hump and approaches with increasing 
temperature the classical value. We therefore expel the 
analysis of d — 1 from the discussions below. For d = 2—6 
the specific heats show humps or peaks, which get sharper 
with increasing d and increasing particle number. How- 
ever, from these smooth curves the orders of the phase 
transition cannot be deduced. 

In Fig. ^| the classification parameters a, 7, t% defined 
above are plotted for two to six dimensions and parti- 
cle numbers up to N = 100. For all values of d the 
parameter a is only a slightly varying function of N and 
approaches very fast an almost constant value. Since a 
is the primary classification parameter from Fig. ^(a) we 
can directly infer that the d = 2 system exhibits a third 
order phase transition (a > 1) while the transition for 
all higher dimensions is of second order (0 < a < 1). 
For N = 50 the dependence of a on d is plotted in 
Fig. |^(a) . Since a decreases rather rapidly with increas- 



5 



(a) 



(b) 



1.5 



1.0 



0.5 



0.0 



0.5 



" o O O O O O o o O o o o o . 



50 100 
N 



150 



0.0 
-0.5 
-1.0 
-1.5 





xxcocccca 


o° 
o 




□ ° 




" o 





50 100 
N 



150 



(c) 



1.00 



0.10 



0.01 



0.00 









Od = 2 






□ d-3 






©d = 4 






Ad = 5 






<d = 6 






10 




100 




N 





FIG. 6: Classification parameters a, 7 and ri for d — 2 — 6 
versus particle numbers N. 



ing d it can be speculated that systems corresponding 
to a large d exhibit a phase transition which is almost 
of first order. As mentioned above for finite systems 
even values a < cannot be excluded by mathemati- 
cal reasons. We note that two-dimensional Bose-gases 
are an interesting and growing field of research. As it 
is well known, the ideal free Bosc-gas in two dimensions 
(d = 1) does not show a phase transition due to ther- 
mal fluctuations which destabilize the condensate p4[ . 
Switching on a confining potential greatly influences the 
properties of the gas, the thermal fluctuations are sup- 
pressed and the gas will show Bose-Einstein condensa- 
tion. Recent experiments have shown that Bose- 
Einstein condensation is possible even though it is called 
a quasi-condensate. In our notion the quasi-condensate 
is just a third order phase transition. Thus, our results 
are in complete agreement with recent experiments and 
earlier theoretical work. An interesting question in this 
respect is whether the order of the transition changes for 
d = 2 in the limit N — > 00. Additional calculation for 
larger N , which are not printed in Fig. |^ indicate that 



a approaches 1 or might even get smaller. Note that 
d = 2 is equivalent to a hypothetical 4-dimensional ideal 
Bose gas or Bosons confined in a 2-dimensional parabolic 
trap. Our results indicate that the order of the phase 
transition sensitively depends on d for values around 2. 
This might be the reason why phase transitions in three 
space dimensions are sometimes classified as second and 
sometimes as third order phase transitions. 

The parameter n is a measure of the finite size of the 
system, i.e. the scaling behavior of t\ as a function of 
A is a measure of how fast a system approaches a true 
n-th order phase transition in the Ehrenfest sense. The 
N dependence of t\ is displayed in Fig. ||(c). The scal- 
ing behavior can be approximated by n ~ N~ s with S 
ranging between 1.06 and 1.12 for d — 2 — 6. 

The d dependence of the classification parameter is 
visualized in Fig. for 50 particles. For this system size 
we found a ~ d~ 4 <^ and t\ 
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FIG. 8: Comparison between calculated zeros of the canon- 
ical partition function for three dimensional trap geometries 
with (a) a continuous single-particle densities of states and 
(b) discrete energy levels for iV = 40. 

The results presented above for continuous single par- 
ticle densities of states ft(E) = E d ~ x are obtained within 
semi- analytical calculations. In order to compare these 
results to systems with a discrete level density we adopt 
as a reference system the 3-dimensional harmonic oscil- 
lator with the partition function given by 



Z{B) = E 



(n + 2)(n + l) 



exp(-B(rc + 3/2)), (12) 



71=0 



with h = lo = &b — 1 • 

Fig. ||(a) displays the zeros of the partition function 
(Eel) for d = 2 and d — 3. Fig. Wh) displays a contour 
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plot of the absolute value of the ground state occupation 
number rj (B) = -±d eo Z(B)/Z{B) with Z given by |[|) 
calculated using an alternative recursion formula . The 
zeros of Z are poles of 770 & n d are indicated by dark spots 
in this figure. 

Analyzing the distribution of zeros consolidates our 
speculation that the order of the phase transition sen- 
sitively depends on d. The distribution of zeros behaves 
like the above calculated values for d — 2 but quanti- 
tatively like d — 3. Since the degeneracy of the three- 
dimensional harmonically confined ideal Bose-gas is a 
second-order polynomial not only the quadratic term has 
to be taken into account. The linear term becomes dom- 
inant for lower temperatures, so for very low tempera- 
tures the best approximation of a continuous one-particle 
density of states is il(E) = E. The parameter a sup- 
ports this statement [Bl, i.e. a resides in a region above 
1. Whereas the parameter 7 behaves like the d — 3 case. 
Finally the parameter t± which is a measure for the dis- 
creteness of the system shows an ~ jy-o.96 dependence 
which is comparable to the one for d — 2. Thus, for small 
systems the phase transition is of third order, it can be 
speculated if it becomes a second order transition in the 
thermodynamic limit. 




400 



FIG. 9: Comparison between three different approaches to 
define a critical temperature for phase transitions in finite 
systems. 



Not only qualitatively but also quantitatively our cal- 
culations are in very good agreement with recent theo- 
retical works [l6|. Comparing the critical temper- 
ature which we defined in Sec. 11 with the usually uti- 
lized ones like the temperature of the peak of the spe- 
cific heat /3(Cv.max) or the grand-canonically calculated 
Tc ~ -ZV 1 / 3 confirms our approach. In Fig. ^ three possi- 
ble definitions of the critical temperature are given which 
all coincide in the thermodynamic limit. All definitions 
show a [3 ~ N~ p dependence with p ranging between 2/5 
and 1/3. 



IV. CONCLUSION 

Starting with the old ideas of Yang and Lee, and Gross- 
mann et al. we have developed a classification scheme for 
phase transitions in finite systems. Based on the ana- 
lytic continuation of the inverse temperature (3 into the 
complex plane we have shown the advantages of this ap- 
proach. The distribution of the so-called Fisher-zeros 
Bk draws enlightening pictures even for small systems 
whereas the usually referred thermodynamic properties 
like the specific heat fail to classify the phase transitions 
properly. The classification scheme presented in this pa- 
per enables us to name the order of the transition in a 
non-ambiguous way. The complex parts Tt of the zeros 
Bk resemble times for which a whole ensemble of identi- 
cal systems under consideration in a heat bath with an 
initial Boltzmann-distribution looses its memory. 

We have applied this to ideal non-interacting Bose- 
gases confined in power-law traps. We have found that 
the order of the phase transition sensitively depends on 
the single particle density of states generated by the con- 
fining potential. The distribution of zeros exactly reveals 
the order of the phase transition in finite systems. 
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